file.loc<-"/Users/ben/evo-dispersal/HZ_Disp/"
fname<-"HZ_disp_mother"

system(paste("R CMD SHLIB ", file.loc, fname, ".c", sep=""))

dyn.load(paste(file.loc, fname, ".so", sep=""))

  # R_n - starting population size
  # R_spX & R_spY - starting habitat matrix size
	# R_lambda - basic number of offspring
	# R_K - carrying capacity
	# R_disp_cost - cost of dispersing .. ranges between 0 & 1
	# R_num_gens - the number of generations past the burn in period... plus also the increase spX to be used
	# rho - the R working environment
.Call('mother',R_num_gens=100, R_n=10000,R_spX=40, R_disp_cost=0.1, R_lambda=5, R_K=100, 
      R_m_rate=0.05, R_maxcost=1, R_beta=2, rho=environment())

inds<-subset(inds, inds[,1]!=-9)
par(mfrow=c(2,1))
plot(inds[,1], inds[,2], col=rgb(red=inds[,4], 0.5, 0.5, 0.5), pch=19)
plot(disp, type="l")

dyn.unload(paste(file.loc, fname, ".so", sep=""))	



### Functions for testing ###

collapse_gaus<-function(x, u) {  
  (exp(-x^2/(2*u^2))/(2*pi*u^2))*(2*pi*x)
}

collapse_gaus_cum<-function(x, u){ # cumulant on the collapse Gaussian
  (u^2-u^2*exp(-x^2/(2*u^2)))/u^2
}

collapse_gaus_rng<-function(n, u){
  r<-runif(n=n, 0, 1)
  x<-u*sqrt(2*log(1/(1-r)))
}

# Function from Barton and gale 1993 - fitness based on Hybrid index
HI_fitness<-function(HI, s, beta){
  W<-1-s*(4*HI*(1-HI))^beta
  W
}